# -*- coding: utf-8 -*-
"""
Created on Wed Aug 24 09:14:13 2022

@author: LE ROUX REMI
"""

# Importation des bibliothèques
import numpy as np
import matplotlib.pyplot as plt

# Données du problème
c = np.array([2e-6, 4e-6, 6e-6, 8e-6, 1e-5])
A = np.array([0.2041, 0.4257, 0.6022, 0.8267, 1.0206])
p_A = 0.025 # précision du spectrophotomètre
u_A = p_A/np.sqrt(3) # incertitude-type sur A (nombre)
u_c = 0.1/100*c # incertitudes-types sur c (liste)

# Régression linéaire
RL = np.polyfit(c, A, 1)
a = RL[0] # coefficient directeur
b = RL[1] # ordonnée à l’origine
print("coefficient directeur : ",a,"L/mol/cm")
print("ordonnée à l’origine : ",b)

plt.plot(c, A, 'go') # données
plt.plot(c, a*c+b,'k--') # droite de régression
plt.xlabel('c/ mol/L'), plt.ylabel('A(631)') # axes
plt.xlim(0, 1.1*np.max(c)) # échelle horizontale
plt.ylim(0,1.1*np.max(A)) # échelle verticale
plt.grid() # fait apparaître la grille
plt.show()

plt.errorbar(c, A, xerr = 2*u_c, yerr = 2*u_A, fmt ='k.')
plt.plot(c, a*c+b,'k--') # droite de régression
plt.xlabel('c/ mol/L'), plt.ylabel('A(631)') # axes
plt.xlim(0, 1.1*np.max(c)) # échelle horizontale
plt.ylim(0,1.1*np.max(A)) # échelle verticale
plt.grid() # fait apparaître la grille
plt.show()

r = A-(a*c+b) # calcul des résidus
plt.errorbar(c, r, yerr = 2*u_A, fmt ='ko')
plt.xlabel('c/ mol/L'), plt.ylabel('résidus')
plt.xlim(0, 1.1*np.max(c)) # échelle horizontale
plt.grid() # fait apparaître la grille
plt.show()

z = (A-(a*c+b))/u_A # calcul des z-scores
plt.plot(c, z,'ko')
plt.xlabel('c/ mol/L'), plt.ylabel('z-score') # axes
plt.xlim(0, 1.1*np.max(c)) # échelle horizontale
plt.ylim(-1.1*np.max(z),1.1*np.max(z)) # échelle verticale
plt.grid() # fait apparaître la grille
plt.show()

# Nombre de simulations
N = 10000

# Simulation des distributions de données
a_sim, b_sim = [], []

for i in range(N):
    A_sim = np.random.normal(A,u_A) # simule A
    c_sim = np.random.normal(c,u_c) # simule c
    RL = np.polyfit(c_sim, A_sim, 1) # regression linéaire
    a_sim.append(RL[0]) # stockage de a simulé
    b_sim.append(RL[1]) # stockage de b simulé

# Calcul des incertitudes-types
u_a = np.std(a_sim, ddof = 1)
u_b = np.std(b_sim, ddof = 1)
print("incertitude-type sur a : u(a) = ",u_a,"L/mol")
print("incertitude-type sur b : u(b) = ",u_b)

#######################################
#######################################

# Importation des bibliothèques
import numpy as np

# Nombre de simulations
N = 1000000

# Données du problème

Vfiole = 100 # mL (volume de la fiole jaugée)
p_fiole = 0.1 # mL (précision de la fiole jaugée)
Vpipette = 10 # mL (volume de la pipette jaugée)
p_pipette = 0.03 # mL (précision de la pipette jaugée)
A631 = 0.7803 # (absorbance de la solution diluée)
p_A631 = 0.025 # (précision de la mesure d'absorbance)
eps631 = 1.017e5 # L/mol/cm (coeff d'absorption molaire)
u_eps631 = 2300 # L/mol/cm (incertitude-type sur le coeff d'absorption molaire)
l = 1 # cm (largeur de la cuve)

# Simulation des distributions de données

Vfiole_sim = Vfiole + np.random.uniform(-p_fiole,p_fiole,N)
Vpipette_sim = Vpipette + np.random.uniform(-p_pipette,p_pipette,N) 
A631_sim = A631 + np.random.uniform(-p_A631,p_A631,N)
eps631_sim = np.random.normal(eps631,u_eps631,N)

c_sim = Vfiole_sim*A631_sim/(Vpipette_sim*eps631_sim*l)

c = np.average(c_sim)
u_c = np.std(c_sim, ddof=1)

print("cmoyen = ",c,"mol/L")
print("u(c) = ",u_c,"mol/L")